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We present the finite amplitude method (FAM), originally proposed in Ref. [I3, for superfluid 
systems. A Hartree-Fock-Bogoliubov code may be transformed into a code of the quasi-particle- 
random-phase approximation (QRPA) with simple modifications. This technique has advantages 
over the conventional QRPA calculations, such as coding feasibility and computational cost. We per- 
form the fully self-consistent linear-response calculation for a spherical neutron-rich nucleus ^^''Sn, 
modifying the HFBRAD code, to demonstrate the accuracy, feasibility, and usefulness of the FAM. 
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I. INTRODUCTION 

Elementary modes of excitation in nuclei provide valu- 
able information about the nuclear structure. The 
random-phase approximation (RPA) based on energy 
density functionals (EDF) is a leading theory applica- 
ble both to low-lying excited states and giant resonances 
P, Although the fully self-consistent treatment of 
the residual (induced) interactions for the realistic energy 
functionals is becoming more and more prevalent |3l-ll3l|. 
the RPA calculations for deformed nuclei are still com- 
putationally demanding. At present, the quasi-particle 
random-phase approximation (QRPA) for deformed su- 
perfluid nuclei are limited only to axially deformed cases 
[l(]| - |l5| . except for Ref. [16] with an approximate treat- 
ment of the pairing interaction. 

Recently, there has been a renewed interest in the so- 
lution of the RPA problem |17-19]. In Ref. the 
finite amplitude method (FAM) was proposed as a feasi- 
ble method for a solution of the RPA equation. The FAM 
allows to calculate all the induced fields as a finite differ- 
ence, employing a computational program of the static 
mean-field Hamiltonian. Recently, the FAM has been 
applied to the electric dipole excitations in nuclei using 
the Skyrme energy functionals fl^- There has been also 
a calculation making use the iterative Arnold! algorithm 
for a solution of the RPA equation [l^. These newly 
developed technologies in conjunction with fast solving 
algorithms for linear systems open the possibility to ex- 
plore systematically the nuclear excitations over the en- 
tire nuclear chart. 

So far, these new techniques (l7l - [l9l | have been devel- 
oped for solutions of the RPA without the pairing cor- 
relations. It is well known, however, that almost all but 
magic nuclei display superfluid features. Therefore, a 
further improvement is highly desirable to make these 
methods applicable to the QRPA equations including 



correlations in the particle-particle and hole-hole chan- 
nels. The purpose of the present paper is to generalize 
the FAM to superfiuid systems, which enables us to per- 
form a QRPA calculation utilizing a static Hartree-Fock- 
Bogoliubov (HFB) code with minor modifications. Our 
final goal would be the construction of a fast computer 
program for the fully self-consistent and triaxially de- 
formed QRPA. This paper is a first step toward the goal, 
to present the basic equations of the FAM for the QRPA 
and show the first result for spherical nuclei. We use the 
spherically symmetric HFB code called HFBRAD [20*1 to 
be converted into the QRPA code. 

This paper is organized as follows: In Sec. [Hi the 
QRPA equation is derived as the small-amplitude limit 
of the time-dependent HFB (TDHFB) equations. In Sec. 
HID we obtain the FAM formulae for the calculation of the 
induced fields. In Sec. IIV[ we summarize all the relevant 
formulae for practical application of the FAM. In Sec. 
|Vl we apply the FAM to the HFBRAD and compare 
the result with that of another self-consistent calculation. 
Sec. IVllis devoted to the conclusions. 



II. SMALL AMPLITUDE LIMIT OF THE 
TDHFB 



In this section, we recapitulate the basic formulation 
of the TDHFB and its small-amplitude limit. In general, 
we will follow the notation in Ref. [l] unless otherwise 
specified. We also use ft = 1 in the following equations. 

We start from the energy functional £[p,K,K*] which 
is a functional of the density matrix and pairing tensor. 
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($|c|cfe|$), Kki = ($|qcfc|$). 
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where |$) is the HFB state. The single-particle Hamilto- 
nian h and the pairing potential A are obtained with a 
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variation of the energy functional with respect to p and 
At*, respectively. 
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The Bogoliubov quasi-particles, (a^,ap, have a lin- 
ear connection to the bare particles, (cfc,c|,); at — 

X]/c(^'=M'^fe + ^fc/iCfe). Here, the index k indicates the 
adopted basis such as the harmonic oscillator states or 
the coordinate space. The quasi-particles are chosen 
so as to diagonalize the HFB Hamiltonian 1]. 
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Here, the normal ordering is assumed. 

In a similar manner, the time-dependent quasi- 
particles flt (i) are characterized by the time-dependent 

wave functions {U{t),V(t)) by at(i) = Y.k{Uk^{t)cl + 
Vkf_i{t)ck}- The time evolution of the quasi-particles un- 
der a one-body external perturbation Fit) are deter- 
mined by the following TDHFB equation. 



^d^^^H{t)+F{tla,{t)], 
where the TDHFB Hamiltonian is given by 
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Here and hereafter, the constant shift is neglected, since 
it does not play any role in the TDHFB equation (U). 
h{t) and A{t) become time-dependent, since they de- 
pend on the densities, p{t) — V*{t)V'^{t) and K{t) — 
V*{t)U'^{t) = -U{t)V^{t), which are time-dependent. 
Note that the static quasi-particles correspond to a quasi- 
static solution of Eq. dH), af,{t) = a^e^-^f'*, with = 0. 

Let us assume that the nucleus is under a weak external 
field of a given frequency uj. 



F{t) = 7?{F(a;)e-*"* + FtHe*"*}, 



-J2fII{cu)b^,, 



(6) 



(7) 



where ^t = at at and i?^^ = at a,y. A small real param- 
eter r] is introduced for the linearization. In the small- 
amplitude limit, the second term (B-part) in Eq. ([7]) 
can be omitted, because it doesn't contribute in the lin- 
ear approximation. The Bogoliubov transformation of 



the external fields {F^^{u!) and F^^{LLi)) is given in Ap- 
pendix IA2I 

The external perturbation F(t) induces a density oscil- 
lation around the ground state with the same frequency 
Lo. The density oscillation, then, produces the induced 
fields in the single-particle Hamiltonian, h{t) = ho + dh{t) 
and in the pair potential, A(i) = A -I- dA{t). Thus, the 
Hamiltonian, Eq. ([S]), is decomposed into the static and 
oscillating parts; H{t) = Hq + 5H{t). 



5H{t) = 'n{SH{uj)e-"^' + SH^{uj)e"^*} 



(8) 



SHiu;) = i^{^i/20H4, + <5i?°2(^)A^^}. (9) 
/it/ 

Here, the i?-part is again neglected in Eq. See Ap- 
pendix [XT] for the derivation of SH{lu). Explicit expres- 
sions for 6H-^^{u}) and 5H^{ijj) are found in Eqs. (jA8|) 
and (|A9p . respectively. 

The time-dependent quasi-particle operators are de- 
composed in a similar manner: 



(10) 



where 5a^(t) can be expanded in the quasi-particle cre- 
ation operators: 

5a^(i)=r;^at (x,,,He— * + i;;He-*) . (11) 



It should be noted that 5a^ can be expanded only in 
terms of the creation operators, because the annihila- 
tion operators in the right-hand side of Eq. (jlip simply 
represent the transformation among themselves, a^(t) = 
X^i/ C^i,{t)av, and do not affect p and k. The amplitudes, 
X and y, must be anti-symmetric to satisfy the fermionic 
commutation relation, {a^{t),a„{t)} — 0. Keeping only 
linear terms in 77, Eq. ^ becomes 



. dSa^{t) 

' dt 



- E^,5a^{t) + [Ho,5a^{t)] + [5H{t) + F{t),a, 



(12) 

Substituting Eqs. (|6l)-(in|) into Eq. ((H]), we obtain the 
linear-response equations: 

{E, + E,+u)Y,A^)+5H%{u) ^F^liu) ' ^'^^ 

In Eq. (jl3p . setting the frequency complex, u u+i^ /2, 
we can introduce a smearing with a width 7. 

Expanding 5H^^{ll!) and SH^'^{uj) in terms of the for- 
ward and backward amplitudes, X and Y, we obtain a 
familiar expression of the equation [1]: 
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This matrix formulation requires us to 
QRPA matrix elements, A^^^fj^i^> and B^ 



fF^"{u)' 
\F''^ij) 

(14) 

calculate the 
This is 



a tedious task and their dimension, which is equal to the 
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number of two-quasi-particle excitations, becomes huge 
especially for deformed nuclei. Instead, in the FAM jl7l |. 
we keep the form of Eq. (fT3|) and calculate the induced 
fields SH^^iuj) and SH'^'^{uj) using the numerical differ- 
entiation. We explain this trick in the next section. 



III. FINITE AMPLITUDE METHOD FOR THE 
INDUCED FIELDS 

The expressions for SH^° and 6H°^ in Eq. are 
given in Eqs. (|A8|) and (|A9p . respectively. Thus, we 
need to calculate Sh{uj) and 6A'^^\uj) for given X and 
Y. We perform this calculation following the spirit of the 
FAM [17]. 

From Eqs. (|10l) and pT|) . we obtain the time- 
dependent quasi-particle wave functions: 



(15) 



where 



Wfe^W = {U + T]{V*X*e'^' + V*Ye-'^')}^ , (16) 



kfi 

Vk^it) = {V + rj{U*X*e'^' + U*Ye~''^')}^^. (17) 

First, let us discuss how to obtain Sh{Lu). The time- 
dependent single-particle Hamiltonian h{t) depends on 
the densities which are determined by the wave func- 
tions {U{t),V{t)). Therefore, h{t) can be regarded as a 
functional of wave functions as 

h [U*{t), V*{t); U{t), V{t)] = h [W{t), V*{t);U{t), V{t)] . 

(18) 

Here, it should be noted that the phase factors, e*^f* 
in Eq. ([T5|) . do not play a role. This is because /i is a 
functional of densities, p, k, and k*, which are given by 
products of one of {U, V) and one of the complex con- 
jugate iU*,V*), such as p = V*V^ and k = V*U^. 
Therefore, the time-dependent phases in Eq. (jlSp are 
always canceled, thus can be omitted. 

Now, we take the small-amplitude limit, keeping only 
the linear order in rj. 

h{t) = h[U*{t),V*{t)-U{t),V{t)] 

= h[U*,V*;U,V]+'n{6h{uj)e-"^^ +\i.c.} .(19) 

Here, 5h{oj) can be obtained using Eqs. ([T6| and ([T7|. 
expanding up to the first order in 77 and collecting terms 
proportional to e"*'^*, as 



dh 

Sh{uj) = ^-VX 
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dh_ 



V*Y 



dh 

dV* 
dh 

^ dV 



■ UX 

■ U*Y. 



(20) 



The calculation of the derivatives, such as dhki/dU^,^, is 
a tedious task and requires a large memory capacity for 
their storage in the computation. In the FAM, we avoid 



this explicit expansion, instead write the same quantity 
as follows: 

^ h[U;,V-U,,V,]-h[U*,V*;U,V] ^ 
V 

(21) 

where U* , V* , Urj , and T/^ are given by 

U: = U*+ T^VX, V* = V*+ r]UX, 



U + rjV*Y, Vr^ = V + r]U*Y. 



(22) 



This is the FAM formula for the calculation of dhiuj). 
All we need in the computer program is a subroutine to 
calculate the single-particle Hamiltonian as a function of 
the wave functions, h[tj* , V*] U, V]. 

For the pair potential, basically, the same arguments 
lead to the FAM formulae for (JA*^^) . The time-dependent 
pair potential A(f) can be written as 

A(i) = A[W{t),V*{t)-U{t)y{t)] 

+ri {(5A(+) (a;)e-''"* + .JA^^^ (w)e''^* } . (23) 

Here, ,5A(+) and JA^^) are independent, since A(i) is 
non-Hermitian in general. (5A(+' can be written in the 
same form as Eq. (PT|) . 



<5A(+)(a;) 



A [C7*, 17*; C/„ K,] - A [[/*, V*-U, V] 



(24) 



where U*, V*, Ur^, and F,, are given by Eq. (f^ . 

The expression for (5A(~) is also obtained from Eq. 
((23)) . collecting terms proportional to e*"*. It is given by 
the same expression as Eq. 



A [u;,v;- [/^, v^] - A [u\v*- u, v] 



(25) 



However, {U*,V*]Uj^,Vn) here are different from Eq. 
([22]) and given by 



u; 



U* + TjVY* 
U + r]V*X* 



; V* + r)UY* 
V + T]U*X* 



(26) 



The essential trick of the FAM is to calculate the in- 
duced fields, (5/i(w) and (5A(='=\ according to Eqs. (PT|) . 
([M]) . and (OSl) with a small but finite parameter 77. Of 
course, the 77^ and higher-order terms bring some nu- 
merical errors, but they are negligible. Therefore, for 
given X and Y , we are able to calculate these induced 
fields, by using the static HFB code with some minor 
modifications. bR'^^(uj) and 6H"^{uj) of Eq. ^ in the 
quasi-particle basis can be calculated with Eqs. (jASp and 
(|A9[) . respectively. Then, we may solve the QRPA linear- 
response equation ([T^ to obtain the self-consistent am- 
plitudes, X and Y, utilizing an iterative algorithm (See 
Sec. lYl). 
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Induced fields in terms of densities 



Eq. (1331), as follows: 



Although the basic formulae of the FAM has been pro- 
vided in Sec. IIIII we may need to modify them in the 
practical implementation of the FAM. For instance, some 
HFB codes, such as HFBRAD, contain subroutines to 
calculate mean fields as functions of densities, not of wave 
functions. In this subsection, we rewrite Eqs. (|2T|) . (l24t . 
and (|25l) in terms of densities. 

The density Sp{t) is written up to linear order in 77 as 



p{t)^V*{t)V^{t) 

= Po + v {Sp{i^)e-'^' + h.c.) 

where 

Sp{uj) = UXV'^ + V*Y^Ul 
This can be written in the FAM form: 

Pn - Po 



(27) 



(28) 



V 



6p{uj) 



where V* and 1^ are given in Eq. (j22l) . 

The pair tensor K{t), which is non-Hermitian, can be 
expressed in a similar manner. 



v*v'^ - v*v'^ 

+ 0{rj'), (29) 



K{t) = V*{t)U^{t) 

Here, k*^^) can be given in the explicit form as 



(30) 



(31) 
(32) 



and in the FAM form as 



+ 0{r,% (33) 



where V* and Ur, are given in Eq. ([22]) for while 

they are given by Eq. (PBl) for 

Now, let us present how to obtain the induced fields 
in terms of the densities. In general, h{t) and /S.{t) may 
depend on p, k, and k*. 

hit) = h [pit), Kit),K*it)] , Ait) = A [pit), Kit),K*it)] . 

(34) 

In order to obtain the induced fields, all we need to do is 
to replace p by defined in Eqs. (P^ . and n by k^^-* in 



Shiio) 



(5A(+)(w) = 



dA'^-\io) = 



(+) (-)- 

Pri, Kri , Kri 



h [p, K, K* 



(+) (-)» 



A[p,K,K*] 



(35) 



(36) 



A 



(-) (+)» 

Pr]^ ^'q 7 ^r] 



A [p, K, K* 



V 



(37) 



where the terms of the second and higher orders in 77 are 
neglected. 



IV. SUMMARY OF THE FINITE AMPLITUDE 
METHOD 

Here we provide a summary of the FAM for the 
QRPA linear-response calculation for a prompt applica- 
tion. Later, we discuss applications of the FAM to the 
Skyrme functionals, however, the FAM formulated in this 
and previous sections is applicable to any kind of energy 
density functional (mean-field) models. 



A. Numerical procedure 

The aim is to solve the linear-response equation 
for a given external field F. In order to obtain the for- 
ward and backward amplitudes, X and Y, we resort to 
an iterative algorithm. Namely, we start from the ini- 
tial guess for (X, F) = iX'-°\Y'-°^) = x<^°\ and calcu- 
late Shiui) and (5A'^^^(ti;) according to the formulae, (j2ip . 
([24| . and ([25]). Then, they are converted into SH'^'^iuj) 
and (5i7"^(aj), using Eqs. (|A8p and (IA9p . respectively. In 
this way, we can evaluate the left and right hand sides of 
Eq. ^ for a given iX,Y). 

Since Eq. (fT3| is equivalent to Eq. (|T4)) . it is a lin- 
ear algebraic equation for the vector x = iX,Y), in the 
form of Ax — b. Many different algorithms are avail- 
able for the solution of linear systems. In this paper, we 
resort to a procedure based on Krylov spaces called gen- 
eralized conjugate residual (GCR) method [21]. Within 
these kinds of methods, a succession of approximate so- 
lutions (if^^^ , x*^^-' , x^^^ , • ■ ■ ) converging to the exact one is 
obtained by the iteration. The GCR algorithm consists 
in a series of steps each containing the operation of the 
matrix A on a given vector, and sums and scalar products 
of two vectors. For the given x = iX, Y), Ax is equal to 
the left hand side of Eq. (fT5)) . Therefore, the quantity 
Ax can be calculated without the explicit knowledge of 
the QRPA matrix itself. 

Here, we summarize the formulae. The linear response 
equation is given by Ax — b, where 



_ (Xfj_^\ r (F]2 
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and 



where 



6HII{lo) = U^ShV* - V^SA^-^*V* 
+ uUa'-+^u* -V^'Sh^U*, 

- V^SA(+^V + U^6h^V. 



Denoting h and A collectively asH = {h, A), the induced 
fields SH are calculated by the FAM formulae, 

^ n[u;,v;;U„v,]~H[u*,v*;u,v] 

where (t/*, V*] U,^,V,j) are given by 

U* = U* + TjVY*, V* = V* + r]UY*, 
= U + r]V*X*, Vri = V + riU* X* . 

for the calculation of Sh{uj) and (5A'^+^. For SA^^\ they 
are 

U* = U* + rjVX, V* = V* + riUX, 

U.r, = U + T]V*Y, V^ = V + T]U*Y. 

The final result does not depend on the parameter rj, 
as far as it is in a reasonable range. The choice of rj is 
discussed in Sec. El 



B. Calculation of the strength function 

Using the solution {X,Y), we can calculate the 
strength function following the same procedure as Ref. 

^Ei^ ^ V \{n\F\Q)\^S{c. E.,) = --lmS{F;u^). 

flf,) TT 



n>0 



(39) 

Here, S{F;uj) is obtained from the solution {X,Y). For 
the operator in the form of Eq. (lAlOp , we may calculate 
S{F;u}) as 

SiF-io) = ti-{pSp{io)}, (40) 
For the operator in the form of Eq. (IA13P , we have 

S{F;u;) = tr \^g'' 6 k^+'^ (uj) + g"< 6 k'--^ (uj)^ . (41) 



For both cases, in the two-quasi-particle basis, Eqs. (|40 
and (|4ip can be written in the unified expression. 



S{F,cj) = i ^ {F'^°*X,Au;) + V(^)} , (42) 



where F^° and are given by Eqs. (jAlip and (|A12p 
for the former case, and by Eqs. (|A14I) and (lAlsp for the 
latter. 



V. APPLICATION OF THE FAM TO THE 
HFBRAD 

In order to assess the validity of the FAM, we install the 
FAM in the HFBRAD code It has to be noted that 
the formalism of the HFBRAD is slightly different from 
the one used in this paper which follows the notations in 
Ref. In particular, the wave functions (iy9i^, (p2/j), the 
pairing tensor p, and the pair potential h are defined in 
a different manner; (pi^{k) = Uk^,, </52^(fc) = V^^, pki = 

s-iid ^ki — ^kiJ where k is the time-reversal state of 
k. A detailed discussion on the difference among the two 
notations can be found in Ref. [25| . 

The HFBRAD ^] is a well known code which solves 
the HFB in the radial coordinate space assuming the 
spherical symmetry. It has been designed to provide fast 
and reliable solutions for the ground state of spherical 
even-even nuclei. For these nuclei, the time-odd densities 
are identically zero and thus they have not been imple- 
mented in the code. In order to render the QRPA fully 
self-consistent, we have to add the time-odd terms in the 
calculation of the induced fields. This task can be sim- 
plified for a case of the presence of spherical and space- 
inversion symmetry, such as in the case of monopole ex- 
citations. For this case, the only time-odd terms with 
non-zero contribution are those due to the current den- 
sity [22], moreover the only non- vanishing component of 
the current density is radial. 

We calculate the strength function of the isoscalar 
monopole for a neutron-rich nucleus, ^'''^Sn. To check 
the self-consistency by looking at the spurious compo- 
nent, we also calculate the strength of the nucleon num- 
ber operator. Both operators are given by the form of 
Eq. (jAlOp with fki = {k\r'^\l) for the isoscalar monopole 
operator and fki = Ski for the number operator. 

In order to obtain the strength function, first, we have 
to solve the HFB equations to construct the ground-state 
wave functions {U,V). It is accomplished by using the 
HFBRAD code. The parameters of the present calcu- 
lation are adjusted to the values used by Terasaki and 
co-workers in Ref. Q ; The box size is Rbox = 20 fm, the 
quasi-particle energy cutoff is E^^ = 200 MeV, the max- 
imum angular momenta of the quasi-particle states are 
Jmax = 21/2 for neutrons, and j^^x = 15/2 for protons. 
We use the Skyrme functional with the SkM* parameter 
set [2^ in the ph-channel and a delta interaction of the 
volume type with the strength Vq = —77.5 MeV fm'^ for 
the pp- and hh-channels. 

The next step is solving the linear-response equation 
for a given external field of the frequency u. At first, we 
build the induced fields, 5h{uj) and 5A'^^\ijj), starting 
from a guess choice of the QRPA amplitudes (A("), F^o'), 
according to Eq. ((38)) . In the present calculation, we 
choose either AT*^*^) = y^") = or the values of A and Y at 
the previous energy cj calculated. We resort to the itera- 
tive algorithm of the GCR method to solve the equation 
([T3|. We include all the two-quasi-particle states {pv) 
within the HFB model space defined above {E^i^^) < 200 
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TABLE I: Convergence properties of the calculation. The 
obtained accuracy e = pf-fe||Vl|6||^ and the number of GCR 
iteration Nuci to reach e < IQ-" are shown for different values 
of 77. The initial vector is chosen as f(") = {X'-°'> ,Y^°^) = (0,0) 
and the maximum number of iterations is set at Mtcr ~ 1, 000. 



MeV). The two-quasi-particle space amounts to 12,632 
states for J'^ = 0+. Note that this number becomes 
much larger if we treat deformed systems. We set the 
accuracy of the convergence to be e < 10"^, where 
e = \\Ax — The number of iterations needed 

depends on u; at low energies, about 50-60 iterations 
are enough to reach the convergence, while, close to the 
central peak at 12 MeV, more than 300 iterations are 
needed. 

We studied the convergence quality of the solutions 
as a function of the parameter rj used for the numeri- 
cal derivative. This is shown in Table HI If 77 is too big 
{ij > 10"'') the derivative of the FAM becomes inaccu- 
rate and the linearity of the procedure is partially broken. 
The residue ||j4a'— &|| reaches a plateau where increasing 
the number of iterations cannot improve it anymore. For 
10"^ < 77 < IQ-^, the calculations converge well and the 
resulting strength function is stable. If rj becomes smaller 
than 10"^", the numerical precision limits are reached 
and the GCR procedure can no longer obtain the required 
precision. Therefore, we may conclude that the parame- 
ter ?7 in the range of 10"^ < ?7 ^ 10-^ is appropriate to 
obtain the induced fields accurately. Although the con- 
stant value 77 = 10"^ is adopted in this paper, we may 
use a more sophisticated choice, such as the w-dependent 
77 values [13, 

We report the strength function of the isoscalar 
monopole mode. To smear the strengths at discrete 
eigenenergies, we add an imaginary term to the energy: 
uj ^ uj + ^7/2, where 7 — 1.0 MeV. This procedure is 
almost equivalent to smearing the strength function with 
a Lorentzian function with a width equal to 7. The cal- 
culated energy- weighted strengths are summed up to 300 
MeV and we found that they exhaust 99.6 % of the the- 
oretical sum-rule value given by ^A{r'^). 

In Fig. [U we compare our results (solid red curve) 
with the one in Ref. [y| (dashed green curve). The self- 
consistent result obtained by Terasaki et al. [6] also em- 
ploys the HFB solutions calculated with the HFBRAD. 
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FIG. 1: (Color onhne) Calculated transition strength of the 
isoscalar monopole 0^ excitations in ^^*Sn (solid red curve), 
compared with the result in |6[ with the cutoff (iii) (green 
dashed curve). The transition strength associated to the 
number operator, magnified by a factor of 10,000, in units 
of MeV-^ is shown by the blue dotted curve. See text for 
details. 



However, in Ref. [6|, the QRPA matrix is calculated 
in the canonical-basis representation and an additional 
truncation of the two-quasi-particle space is introduced 
for the construction of the QRPA matrix . In contrast, 
we introduce no additional truncation for our FAM cal- 
culation. We compare our results with the one of the 
cutoff (iii) in Ref. Q which takes into account the high- 
est number of states for the construction of the QRPA 
matrix; all the proton quasi-particles up to 200 MeV and 
the neutron canonical levels with occupancy > 10"^^. 

In the first two peaks at ~ 5 and 8.5 MeV, the two 
curves are almost perfectly overlapping. The peaks be- 
tween 11 MeV and 18 MeV occur at the same energy 
for the two calculations while their height is slightly dif- 
ferent. The bump close to zero energy resulting in our 
calculations has to be attributed to the presence of a spu- 
rious mode. To check the position of the spurious mode 
related to the pairing rotation of the neutrons, we in- 
cluded in Fig. [T]the transition strength associated to the 
number operator, by the blue dashed line. The spurious 
mode is well localized close to zero energy. 

The present result demonstrates the accuracy and use- 
fulness of the FAM for the superfluid systems. Even if the 
two codes include some differences in the truncation of 
the two-quasi-particle space, the similarity of the results 
is very satisfying 



VI. CONCLUSIONS 

The finite amplitude method for the QRPA has been 
presented. The basic idea is identical to the original FAM 
[17|, that we resort to a numerical differentiation to calcu- 
late the induced fields and then solve the linear-response 
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equation with an iterative algorithm such as the GCR. 
With the FAM, a HFB code with simple modifications 
can be turned into a QRPA code. Especially, it is very 
easy to construct the QRPA code which has the same 
symmetry of the parent HFB one whose subroutines are 
used to perform the numerical derivative. All the terms 
present in the TDHFB calculation, including the time- 
odd mean fields, should be taken into account to con- 
struct fully self-consistent codes. This requires us some 
effort to update the original HFB code. Still, the neces- 
sary task for coding the FAM is much less than that for 
the explicit calculation of the QRPA matrix elements for 
realistic energy functionals. In addition, it does not re- 
quire a large memory capacity, since we do not construct 
the QRPA matrix. We have built a fully self-consistent 
QRPA code using the HFBRAD [20|]. The iterative al- 
gorithm, for which we adopted the GCR method in this 
paper, may be replaced by a better algorithm in future. 
The resulting strength functions of the isoscalar 0'*' mode 
of ^^^Sn show high similarity with the fully self-consistent 
calculations in Ref. @. Thus, this paper showed the 
first application of the FAM for superfiuid systems and 
demonstrated the usefulness of the FAM for the construc- 
tion of the QRPA code by modifying existing HFB codes. 
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Here, Sh(t) and SA{t) are oscillating as 

Sh{t) = 77(,5/i(a;)e^'"* + ,5/i^(a;)e*"*) , (A2) 
SAit) = 7] (5A(+)(a;)e-''"* + SA^-\uj)e''^'^ . (A3) 

Note that SA'^^^uj) are anti-symmetric but Sh{uj) is not 
necessarily Hermitian. The induced Hamiltonian, Eq. 
(|A1[) . is now expressed in the form of Eq. ([SJ with 6H{uj) 
given by 




Hereafter, Sh{uj) and 5A^-^\u}) are denoted by Sh and 
SA^^\ for simplicity. 

Since the Bogoliubov transformation can be written in 
terms of the unitary matrix W [l| as follows: 

-'(;)• 

we may rewrite Eq. (IA4[) in the quasi-particle basis: 

^ \ («' «) - (-ii.. 1") - (:) ^ 

(A6) 

This transformation should provide 5H^° and dH'^^ in 
Eq. dH). 

_^SH^YJ \-SA(-> -5h^J ■ 

(A7) 

We write here their explicit expression: 

SH^liij) = (^U^ShV* -V^6A'^->V* 

+U''SA^+'>U* - V^'Sh'^U*) , (A8) 

6H"^l{oj) = (^-V^ShU + U^6A(->U 

-V^6A(+^V + U^Sh^v) . (A9) 



Appendix A: Bogoliubov transformation of 
one-body fields 

1. Induced fields SH 

The TDHFB Hamiltonian is given by Eq. ©. We 
consider the small-amplitude limit, H{t) = Hq + SH{t), 
where Hq is the HFB Hamiltonian of Eq. 1^ and 



SH{t) 



Sh{t) 
-SA*{t) 



6A{t) 
-5h*(t) 




(Al) 



2. External field F 

The one-body field in general can be written in a form 
of Eq. ([7]) in terms of the quasi-particle operators, ne- 
glecting a constant. Suppose that F{lij) in Eq. ([5]) has a 
form 



(AlO) 



where the difference of a constant shift is neglected. Here, 
the matrix fki is a general complex matrix, since F{u)) is 
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non-Hermitian in general. The Bogoliubov transforma- the same calculation provides F'^^ and F'^^ by 
tion as in Eq. (IX7|) . then, leads to F^" and i^"^ in Eq. 



Fll = {U^fV*-V^fU*)^^, (All) 

FHI = (U^fV-V^fU)^^. (A12) ^ :^.r, 



In case that -F(w) has a form of pairing-type 



(.9fe'44 + g'u^ici) , (A13) 



F'" = {U^gU*-V^g'V*)^^^ (A14) 
^^^^ ^ (U^g'U-V^gV)^^. (A15) 
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